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I.  INTRODUCTION 
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A  new  approach  is  described  for  combining  range  and 
Doppler  data  from  multiple  radar  platforms  to  perform 
multi-target  detection  and  tracking.  In  particular,  azimuthal 
measurements  are  assumed  to  be  either  coarse  or  unavailable, 
so  that  multiple  sensors  are  required  to  triangulate  target  tracks 
using  range  and  Doppler  measurements  only.  Increasing  the 
number  of  sensors  can  cause  data  association  by  conventional 
means  to  become  impractical  due  to  combinatorial  complexity, 
i.e.,  an  exponential  increase  in  the  number  of  mappings 
between  signatures  and  target  models.  When  the  azimuthal 
resolution  is  coarse,  this  problem  will  be  exacerbated  by  the 
resulting  overlap  between  signatures  from  multiple  targets  and 
clutter.  In  the  new  approach,  the  data  association  is  performed 
probabilistically,  using  a  variation  of  expectation-maximization 
(EM).  Combinatorial  complexity  is  avoided  by  performing 
an  efficient  optimization  in  the  space  of  all  target  tracks  and 
mappings  between  tracks  and  data.  The  full,  multi-sensor, 
version  of  the  algorithm  is  tested  on  simulated  data.  The 
results  demonstrate  that  accurate  tracks  can  be  estimated  by 
exploiting  spatial  diversity  in  the  sensor  locations.  Also,  as  a 
proof-of-concept,  a  simplified,  single-sensor  range-only  version 
of  the  algorithm  is  tested  on  experimental  radar  data  acquired 
with  a  stretch  radar  receiver.  These  results  are  promising,  and 
demonstrate  robustness  in  the  presence  of  nonhomogeneous 
clutter. 
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We  present  a  new  approach  for  multi-target 
detection  and  tracking,  in  which  information  from 
multiple,  spatially  diverse,  radar  sensors  is  combined 
to  improve  track  reliability  and  accuracy.  In  particular, 
we  treat  the  difficult  case  in  which  azimuthal 
measurements  are  either  coarse  or  unavailable,  so 
that  multiple  sensors  are  required  to  triangulate 
target  tracks  using  range  and  Doppler  measurements 
only.  However,  data  association  in  this  problem  is 
extremely  complex  for  several  reasons.  First,  the 
coarse  azimuthal  resolution  can  result  in  significant 
overlap  between  the  signatures  of  the  multiple 
targets  and  clutter.  Second,  increasing  the  number 
of  sensor  platforms  leads  to  an  exponential  increase 
in  the  number  of  mappings  between  signatures  and 
hypothesized  targets.  Thus,  it  becomes  impractical 
to  sort  out  the  associations  using  combinatorial 
approaches. 

Multiple  hypothesis  tracking  (MHT)  [1,  2],  a 
benchmark  multi-target  tracking  algorithm,  performs 
data  association  using  an  exhaustive  evaluation  of 
all  mappings  between  targets  and  data  samples,  and 
is  therefore  subject  to  a  combinatorial  explosion 
as  the  amount  of  data  and  the  number  of  sensors 
increase.  Pruning  or  gating  are  typically  used  to 
alleviate  the  computational  burden  by  eliminating 
the  less  likely  hypotheses,  however  valid  hypotheses 
may  also  be  eliminated  in  the  process.  Therefore, 
data  having  low  signal-to-clutter  ratio  (S/C)  may 
be  discarded,  which  can  lead  to  missed  detections. 

An  alternative  to  MHT  is  joint  probabilistic  data 
association  (JPDA)  [3-5]  which  is  more  efficient 
than  MHT  because  one  only  needs  to  evaluate  the 
association  probabilities  separately  at  each  time  step. 
Whereas  more  simplistic  single  scan  methods  (such 
as  nearest  neighbor  approaches)  consider  only  the 
single  observation  closest  to  the  predicted  state,  JPDA 
is  more  robust  since  the  state  of  each  track  is  updated 
using  a  weighted  average  of  all  measurements  falling 
within  its  validation  region  at  the  current  time  step. 
However,  since  it  is  a  single  scan  process,  not  all 
possible  data-to-track  mappings  are  considered.  A 
drawback  of  JPDA  is  that  while  it  is  appropriate  for 
track  maintenance,  it  lacks  an  explicit  mechanism  for 
track  initiation  [6].  Both  MHT  and  JPDA  have  been 
adapted  for  multi-sensor  scenarios  [5]. 

A  data  association  approach  based  upon  linear 
programming  has  been  proposed  which,  like  JPDA, 
updates  the  track  states  using  weighted  averages  of 
measurements  [6].  Results  from  computer  simulations 
indicate  significantly  lower  computational  complexity 
than  JPDA,  as  well  as  improved  accuracy.  Moreover, 
this  method  provides  an  explicit  mechanism  for  track 
initiation. 

Recently,  sequential  Monte  Carlo  methods,  a.k.a. 
“particle  filters,”  have  been  adapted  for  multi-target 
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tracking  problems  [7-9].  Whereas  Kalman  filters 
are  used  in  traditional  JPDA  methods  for  updating 
track  states,  particle  filters  are  more  appropriate  for 
situations  involving  nonlinear  state  and  measurement 
equations  and  non-Gaussian  noise.  Data  association 
for  particle  filter  methods  can  be  performed  in  various 
ways,  as  described  in  [8]. 

Both  MHT  and  JPDA  assume — often 
correctly — that  a  target  can  generate  at  most  one 
measurement  per  scan.  However,  if  this  constraint 
is  relaxed,  data  association  can  be  formulated 
using  a  continuous  optimization  procedure,  notably 
expectation-maximization  (EM)  [10-12],  rather  than 
combinatorics.  Thus,  an  efficient  “hill  climbing” 
optimization  is  performed  in  the  space  of  all  model 
parameters  and  all  possible  mappings  between  data 
samples  and  hypothesized  targets.  An  important 
advantage  of  this  formulation  is  that  computational 
complexity  scales  only  linearly  with  the  number 
of  targets  and  sensors,  whereas  combinatorial 
approaches  such  as  MHT  scale  exponentially.  Thus, 
trackers  based  upon  the  EM  formulation  would,  in 
principle,  be  practical  for  scenarios  with  high  clutter 
and/or  high  densities  of  targets.  Streit,  et  al.  [13-15] 
developed  the  probabilistic  multihypothesis  tracker 
(PMHT)  that  utilizes  EM  to  perform  data  association 
while  simultaneously  estimating  tracks  based  upon 
multiple  scans  of  data.  This  approach  has  also  been 
extended  by  other  authors,  for  example  to  cases 
with  multiple  sensors  [16-18]  and  maneuvering 
targets  [19].  Avitzour  [20]  and  Perlovsky  [21,  22] 
also,  independently,  developed  maximum-likelihood 
procedures  for  multi-target  tracking  which  use 
EM  for  data  association.  Perlovsky  demonstrated 
mathematically  [22,  23],  using  Cramer-Rao  bound 
analysis,  that  the  utilization  of  classification  features 
within  the  tracker  is  equivalent  to  an  improvement 
in  S/C  ratio  in  high-clutter  tracking  environments. 
Subsequently,  he  developed  a  version  of  the  algorithm 
in  which  classification  and  tracking  are  performed 
concurrently  [22,  24-26].  Here,  if  they  are  available, 
classification  features  (e.g.,  radar  cross  section  (RCS), 
length,  etc.)  are  placed  on  an  equal  footing  with 
tracking  features  (e.g.,  range,  Doppler,  bearing,  etc.). 
The  model  is  then  a  mixture  of  different  types  of 
target  and  clutter  components  in  the  combined  space 
of  tracking  and  classification  features.  Like  PMHT, 
Perlovsky’ s  approach  is  a  multi-scan  (i.e.,  “batch”) 
algorithm. 

Perlovsky ’s  approach  differs  from  PMHT  chiefly 
in  the  choice  of  track  model.  In  PMHT  the  motion 
of  each  target  is  modeled  using  a  set  of  discrete-time 
state  transition  equations,  and  therefore  a  Kalman 
smoother  is  used  to  estimate  track  parameters 
in  the  M-step  of  the  EM  iterations  [13-15]  (in 
nonlinear/non-Gaussian  cases  the  target  states  are 
obtained  using  dynamic  programming  rather  than 
Kalman  smoothers).  In  contrast,  Perlovsky’s  approach 


is  flexible  with  respect  to  the  choice  of  track  model, 
but  normally  the  choice  is  to  include  continuous 
polynomial  models  (e.g.,  constant  velocity,  constant 
acceleration,  etc.),  or  piecewise  polynomial  models, 
for  target  trajectories.  The  use  of  polynomial  models 
leads  to  very  simple  parameter  update  formulas  in  the 
M-step,  for  example,  simple  matrix  inversions  which 
are  similar  in  structure  to  polynomial  regression  [22]. 
Of  course,  under  certain  conditions  the  discrete-time 
state  transition  equations  used  in  PMHT  would  be 
equivalent  to  polynomial  models.  Another  difference 
is  in  the  choice  of  optimization  criterion.  Whereas  in 
PMHT  the  goal  is  maximization  of  the  a  posteriori 
probability  (MAP)  [15],  Perlovsky  uses  maximum 
likelihood  estimation  (MLE)  in  the  case  of  sampled 
data,  and  minimization  of  the  cross-entropy  in  the 
case  of  pixelated  data  [22].  Avitzour’ s  procedure 
is  maximum  likelihood,  like  Perlovsky’s,  however 
different  tracking  models  are  specified  due  to 
differences  in  the  particular  applications. 

The  algorithm  described  in  the  present  paper  is 
developed  along  the  lines  of  Perlovsky' s  general 
approach,  and  utilizes  EM  for  data  association. 
However,  a  major  new  aspect  considered  here  is 
the  data  model,  which  is  particularly  complicated, 
and  arises  from  combining  data  from  multiple 
sensors,  where  azimuthal  measurements  are  absent. 
Thus,  the  goal  here  is  to  use  multiple  sensors  to 
triangulate  target  tracks  using  range  and  Doppler 
measurements  only,  while  performing  data  association 
probabilistically.  An  abbreviated  version  of  this  work 
was  published  previously  [30,  31],  although  some 
results  and  most  of  the  mathematical  details  were  not 
included. 

An  important  consideration  for  any  tracking 
algorithm  is  computational  cost.  When  applied  to 
a  benchmark  (single-sensor,  single-target)  tracking 
problem,  it  was  found  that  the  computational 
cost  of  PMHT  is  roughly  the  same  order  of 
magnitude  as  the  cost  of  MHT  and  JPDA  [27],  The 
computational  cost  of  PMHT  versus  other  methods 
has  also  been  evaluated  for  multi-target  tracking 
[28,  29].  Perlovsky’s  approach  would  likely  have 
a  computational  cost  similar  to  PMHT,  since  they 
share  a  similar  structure,  and  are  both  based  upon 
EM.  For  multi-sensor,  multi-target  applications, 

MHT  and  JPDA  would  likely  require  hypothesis 
pruning  in  order  to  avoid  an  exponential  increase 
in  computational  cost.  Since  the  computational  cost 
of  EM-based  approaches  scale  only  linearly  with 
increasing  numbers  of  sensors,  it  is  expected  that  the 
advantages  of  EM-based  approaches  would  become 
increasingly  evident  for  multi-sensor  applications, 
such  as  the  one  considered  in  this  paper. 

The  paper  is  organized  as  follows.  In 
Sections  II-IV  the  mathematical  approach  is 
developed  for  the  full,  multi-sensor,  version  of  the 
algorithm.  The  EM  algorithm  is  a  well-established 


594  IEEE  TRANSACTIONS  ON  AEROSPACE  AND  ELECTRONIC  SYSTEMS  VOL.  45,  NO.  2  APRIL  2009 


Authorized  licensed  use  limited  to:  AFRL.  Downloaded  on  July  13,  2009  at  14:21  from  IEEE  Xplore.  Restrictions  apply. 


mathematical  tool  for  estimating  the  parameters  in 
mixture  models  [10-12],  and  could  certainly  be 
used  as  a  framework  for  deriving  the  solution  to 
our  problem,  as  was  done  for  the  PMHT  [13-15]. 
However,  in  this  paper  a  different  approach  is 
offered,  which  Redner  and  Walker  [12]  refer  to  as 
the  “traditional  general  approach.”  Here,  a  system 
of  likelihood  equations  is  derived,  then  solved  using 
the  partial  derivatives  of  the  optimization  criterion, 
as  described  in  [22],  [32].  This  derivation  does  not 
require  an  understanding  of  EM.  In  addition  to  the 
full,  multi-sensor,  algorithm,  a  simplified  version 
is  also  derived  in  Section  V  which  is  appropriate 
for  performing  range-Doppler-only  tracking  from 
single-sensor  data.  In  this  case,  the  update  formula 
for  track  parameters  is  particularly  simple  and 
efficient — it  consists  of  a  small  matrix  inversion  for 
each  target  component.  Most  of  the  mathematical 
details  are  contained  in  the  appendices.  Appendix 
A  contains  a  derivation  of  the  parameter  estimation 
equations  for  general  mixture  models,  and  also 
tailors  these  equations  to  the  tracking  model  used 
in  this  paper.  Appendix  B  contains  a  convergence 
proof  which,  again,  is  within  the  structure  of  the 
traditional  general  approach  rather  than  depending 
upon  EM.  Finally,  in  Section  VI  results  are  presented 
from  experiments  designed  to  test  the  algorithm. 

The  simplified,  single-sensor  version  is  tested  on 
experimental  radar  data,  while  the  full,  multi-sensor, 
version  is  tested  on  synthetic  data.  The  results, 
although  preliminary,  demonstrate  robustness  in  the 
presence  of  nonhomogeneous  clutter,  and  uncertainty 
in  the  number  of  targets  present.  Section  VII  provides 
a  summary  and  discusses  directions  for  future 
research. 

II.  DESCRIPTION  OF  THE  SENSORS  AND  DATA 

The  sensor  model  for  this  discussion  incorporates 
a  ground-based  radar  antenna  having  poor  azimuthal 
resolution.  In  fact,  we  will  assume  the  extreme  case  in 
which  each  sensor  measures  target  range  and  Doppler 
(range-rate)  only,  but  no  azimuth.  Thus,  the  track  of 
each  target  can  only  be  estimated  by  triangulation 
from  multiple,  spatially  diverse,  sensor  platforms.  In 
this  work,  we  focus  on  the  2-dimensional  tracking 
problem  in  which  targets  are  assumed  to  lie  on  the 
zero-elevation  plane,  while  sensors  have  arbitrary 
3-dimensional  positions. 

Our  method  is  appropriate  for  any  collection  of 
range/Doppler  data,  although  as  a  working  model 
it  is  assumed  the  data  are  acquired  using  a  stretch 
receiver  [33].  Here,  a  suite  of  long-duration  chirps 
is  transmitted  from  a  stationary  ground-based  radar, 
and  the  corresponding  suite  of  received  signals  within 
the  coherent  processing  interval  (CPI)  is  processed 
jointly  to  produce  a  two-dimensional  digitized  image 
in  range/Doppler  coordinates,  referred  to  as  a  “scan 


Fig.  1.  Sample  stack  of  three  pixelated  scan  frames  acquired 
with  a  particular  sensor.  Target  signature  appears  as  high  intensity 
(dark)  blob  centered  around  its  true  range/Doppler  coordinate. 
Since  target  is  in  motion,  signature  position  is  slightly  different  in 
each  scan  frame.  Similar  data  stacks  would  be  acquired  using 
other  sensors,  and  entire  set  of  stacks  from  all  sensors  is 
processed  jointly  to  estimate  target  tracks. 

frame.”  Examples  of  scan  frames  are  shown  in  Figs.  2 
and  3  in  Section  VI,  where  the  target  and  clutter 
signatures  appear  as  blobs  of  energy  centered  upon 
their  true  range  and  Doppler.  The  brightness  (strength) 
of  each  blob  is  proportional  to  the  range  and  RCS 
of  the  reflectors.  The  target  and  ground  clutter  may 
occupy  several  range-Doppler  resolution  cells.  Of 
course,  the  range  resolution  is  inversely  proportional 
to  the  bandwidth  of  the  transmitted  signal,  and 
is  degraded  by  the  associated  pulse  compression 
processing.  Similarly,  the  Doppler  resolution  is 
inversely  proportional  to  the  CPI,  and  is  degraded 
by  sensor  platform  vibration,  system  instability,  and 
processing.  The  spreading  of  the  target  signature 
across  several  range-Doppler  cells  may  be  associated 
with  target  speed  and  acceleration,  target  size,  and 
scintillation  of  target  cross  section.  Ground  clutter 
appears  in  the  scan  frame  image  as  a  ridge  centered  at 
zero  Doppler  and  spread  across  all  range  bins  (again, 
refer  to  Figs.  2  and  3  in  Section  VI).  The  spread  of 
the  ridge  beyond  the  Doppler  resolution  cell  may  be 
caused  by  antenna  scan  modulation,  sensor  platform 
vibration,  as  well  as  internal  motion  of  the  clutter 
(e.g.,  movement  of  leaves  and  branches  in  the  wind). 
Depending  on  its  width,  this  clutter  ridge  may  in  some 
cases  obscure  or  partially  obscure  the  signatures  from 
targets  with  small  radial  velocities.  Finally,  there  will 
be  a  certain  amount  of  background  noise  (receiver 
noise)  which  is  uniformly  distributed  in  range/Doppler 
over  the  image. 

It  is  assumed  the  data  are  collected  as  follows 
(Fig.  1).  At  each  time  tj,  for  j  =  1,2,3,.  ..,7,  the 
received  signals  within  the  CPI  centered  at  tj  are 
processed  jointly  to  yield  a  sampled  or  pixelated 
range/Doppler  image  (scan  frame)  as  described  above. 
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Thus,  for  each  sensor  m  =  1,2,3 we  acquire  a 
stack  of  J  range/Doppler  scan  frames  corresponding 
to  the  different  times  t=.  Note  that,  due  to  the  irregular 
overlap  in  the  fields-of-view  of  the  different  sensors, 
the  signatures  from  some  targets  may  not  appear 
in  the  images  from  all  sensors  at  all  times.  Each 
range/Doppler  scan  frame  image  is  described  using 
the  notation  p0(wymn),  where  p0  is  the  pixel  intensity 
at  the  range/Doppler  coordinate  wjmn  =  (rjmn,djmn), 
which  is  indexed  by  time  index  j  (i.e.,  scan  frame 
number  j),  sensor  m,  and  pixel  n  =  1,2,3 
The  total  number  of  pixels  in  each  scan  frame  is 
N  =  Nr  x  Nd,  where  Nr  and  Nd  are  the  number  of  bins 
in  range  and  Doppler,  respectively.  The  width  of  each 
pixel  with  respect  to  range  and  Doppler  is  Ar  and 
Ad,  respectively.  The  data  from  the  multiple  sensor 
platforms  are  combined  noncoherently  to  estimate 
target  tracks,  as  described  in  the  following  sections. 

III.  MODEL  FOR  THE  DATA 


As  described  in  Section  II,  the  target  signatures 
appear  as  blobs  of  energy  in  the  scan  frame  images 

Pofrjmn)’  where  reca11  that  Wjm„  =  (■ rJmn’djmn )  is 
the  range  and  Doppler  value  of  the  nth  pixel  in 
the  frame  pertaining  to  the  mth  sensor  at  the  j th 
time.  The  signature  corresponding  to  target  k  is 
roughly  centered  at  the  true  range  R-km  and  Doppler 
D  ,  ,  measured  relative  to  sensor  m  at  time  t,. 

We  wish  to  construct  a  mathematical  model  that 
approximates  these  signatures,  and  for  target  k  we 
denote  this  model  as  pk(v/jmn  \  Qk),  where  Qk  is  the 
set  of  model  parameters.  In  this  case  it  makes  sense 
to  use  Gaussian  distributions  to  model  the  target 
signatures  (blobs  of  energy),  so  that  the  signature 
due  to  target  k  for  sensor  m  and  time  t f  is  modeled 
by 


Pk^imn  I  0a)  = 


A,- A, 

2™ 7rkadk 


exp  < 


r  —  R 

jnm  jkm 


^ jnm  D  jkm 


A  model  for  the  data  p0(w/m„)  can  be  developed 
which  depends  upon  the  target  trajectories  as  well 
as  the  sensor  parameters  and  coordinates.  Suppose 
the  coordinates  of  target  k  at  time  t  -  are  given  by 
the  (east,  north)  coordinate  (xk(tj),yk(tj))  =  (xjk,yjk). 
Using  a  constant  acceleration  model,  the  trajectories 
are  described  by 


Xjk=Xk+x'ktj+4tf 

(i) 

ykk  =  yk  +  y’ktj  +  A2- 

(2) 

Here  (xa,ya)  is  the  time-zero  position  of  target  k, 
while  (x'k , y'k )  is  the  time-zero  velocity  and  (xk,yk) 
are  proportional  to  the  x  and  y  components  of  the 
acceleration.  The  targets  are  assumed  to  lie  at  zero 
elevation,  which  simplifies  the  discussion,  as  well  as 
reducing  the  number  of  track  parameters  that  need 
to  be  estimated.  The  sensors  are  allowed  to  have 
arbitrary  3-dimensional  coordinates.  Suppose  there 
are  M  sensor  platforms,  where  sensor  m  is  fixed  at  the 
(east,  north,  elevation)  coordinate  (xm,ym,zm).  Then, 
the  range  Rkm(tj )  =  Rjkm  from  sensor  m  to  target  k  at 
time  tj  is  given  by  the  equation 

Rjkm  =  \j (*m  -  Xjk)2  +  (ym  -  yjk y-  +  4-  (3) 

The  Doppler  (range-rate)  Dkm(tj)  =  [dRkm(t)/dt\t=tj  is 
then 

DknMj)  =  D  jkm 


(5) 

Here,  the  set  of  model  parameters  is  <dk  = 
Wdk^rkA°k,yk,x'k,y'k,xk,yk}  Where,  in  (5),  the 
track  parameters  {xk ,  yk  ,xk,yk,  xk ,  yk  }  are  contained 
implicitly  within  the  quantities  R-km  and  Djkm 
according  to  (3)  and  (4).  Also,  Ar  and  Ad  define 
the  pixel  size  relative  to  the  range  and  Doppler 
coordinates,  as  discussed  in  Section  II.  Note 
that  the  distribution  in  (5)  is  normalized1  so  that 
Yfn=  i  Pk^jmn  I  ©*)  =  1-  The  variance  parameter 
er^.  specifies  the  spread  (width)  of  the  signatures 
with  respect  to  the  range  coordinate  in  the  scan 
frame  image,  and  is  related  to  the  range  resolution 
of  the  sensor  and  the  range  extent  of  the  target. 

The  other  variance  parameter  adk  specifies  the 
width  of  the  signature  in  Doppler,  and  depends 
both  on  the  sensor  resolution  and  the  target  motion, 
including  internal  (micro)  motion.  It  should  be 
noted  that  the  model  could  alternatively  have 
been  formulated  so  that  the  variance  parameters 
depend  upon  time  index  j  or  sensor  m ,  as  well  as 
sensor  k. 

The  interference  also  requires  a  model,  and  we 
approximate  this  using  two  components,  one  for 
the  uniform  background  noise  (receiver  noise),  and 
one  for  the  stationary  clutter.  We  designate  the  two 
interference  indices  as  k  =  (K  —  1),  K  and  the  target 
indices  as  k  =  1,2,3,. . .,(K  —  2).  The  background 
noise  (receiver  noise)  component  k  =  K  —  1,  uniform 
(constant)  in  both  range  and  Doppler,  is  described  by 


.v„ ,  —  a: 


jk 


R 


jkm 


(x'k  +  2 4't  •) 


Pk-\ (w/m«  I  ©X-l)  - 


(6) 


(  ym-yjk 
V  Rjkm 


( v'k  +  Zy'ktj )■ 


(4) 


'Actually,  the  normalization  is  approximate,  and  assumes  the  pixel 
dimensions  (A;.,Arf)  are  small  relative  to  (<rrk,adk). 
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where  N  is  the  number  of  pixels  in  each 
range/Doppler  image,  as  discussed  above.  The 
1  /N  factor  is  necessary  for  normalization,  i.e.,  to 
insure  J2n=iPk(wjmn  I  0*)  =  1-  Note,  here  the  set  of 
parameters  (-)K  ,  is  empty.  The  other  interference 
component  k  =  K,  which  models  the  reflected  energy 
from  stationary  ground  clutter,  is  uniform  (constant)  in 
range,  and  zero-mean  Gaussian  in  Doppler,  i.e.. 


PK&imn  I  Qk)  = 


A, 


V2^adkNr 


exp 


(7) 


Here,  Nr  is  the  number  of  range  bins  in  each  pixel 
array,  as  discussed  in  the  previous  section.  Note, 
here  the  set  of  parameters  is  0A-  =  { rrjjK } .  A  further 
discussion  of  clutter  modeling  is  given  in  [22], 

The  total  model  for  the  image  data  is  the  weighted 
summation  (mixture)  of  individual  target  and  clutter 
components 

K 

Pfrjmn  I  0)  =  Y,EkmPk('/Vjmn  I  0i)-  (8) 

k=  1 


we  will  investigate  how  model  selection  strategies 
described  in  the  literature  can  be  adapted  to  the 
tracking  problem,  for  choosing  the  number  of  target 
components,  pruning  unneeded  components,  and 
adding  components  when  necessary. 

A  useful  result  can  be  obtained  by  summing  both 
sides  of  (8)  over  the  pixel  index  n,  and  using  the 
(previously  stated)  fact  that  £^=  i  Pk(wjm„  I  0*)  =  1  > 
to  obtain 

Y.P^jmn  I  0)  =  YsEkm  (9) 

n  k 

for  each  j  =  1,2 ,...,/.  Here  we  have  assumed  that  if 
a  target  is  in  the  field-of-view  of  a  particular  sensor 
m,  it  will  remain  in  the  field-of-view  for  the  entire 
set  of  scan  frames  j  =  1,2,..., 7.  Thus,  the  sum  of 
the  relative  model  weights  is  independent  of  j.  It 
is  sensible  to  place  a  constraint  so  that,  for  each 
scan  frame  j  and  sensor  m,  the  total  energy  in  the 
model  equals  the  total  energy  in  the  measured  data, 
i.e.,  £„p(wjm„  I  0)  =  £nPo(wym«)-  Combining  this 
constraint  with  (9),  we  obtain  the  following  constraint 
on  the  mixture  weight  parameters,  valid  for  every  scan 
frame  j  and  sensor  m: 


Here,  Ekm  is  the  relative  weighting  for  each  of  the 
model  components,  which  we  refer  to  as  the  mixture 
weight.  For  target  components,  the  mixture  weight  is 
proportional  to  the  target’s  RCS  and  range  as  viewed 
by  sensor  m.  Since  RCS  can  be  strongly  dependent 
upon  aspect  angle,  we  allow  Ekm  to  depend  upon  both 
target  k  and  sensor  m.  This  dependence  upon  sensor 
ni  also  allows  for  the  fact  that  the  fields-of-view  for 
the  different  sensors  overlap  in  an  irregular  manner,  so 
that  a  signature  from  a  certain  target  may  not  appear 
in  the  data  from  all  sensors.  In  (8),  the  total  set  of 
model  parameters  is  denoted  by  0  =  { Ekm ,  (-)k } ,  where 
k  =  1,2 and  m  -  1,2 ,...,M. 

The  well-known  problem  of  estimating  the  best 
number  of  components  K  in  mixture  models  is 
discussed  in  the  literature.  For  example,  mathematical 
strategies  are  described  in  [34],  [35],  [36]  and  the 
references  therein.  In  our  application,  this  so-called 
“model  selection”  problem  stems  from  the  fact  that  the 
true  number  of  targets  is  not  known  in  most  contexts, 
and  therefore  we  must  make  an  educated  guess  for 
the  appropriate  number  of  model  components  K. 

If  there  are  more  model  components  than  targets, 
unneeded  components  can  presumably  be  eliminated 
automatically  as  their  mixture  weights  Ekm  adapt 
to  small  values.  A  tracking  example  is  presented  in 
Section  VI  where  this  type  of  process  is  illustrated. 

In  cases  where  the  true  number  of  targets  exceeds 
the  number  of  model  components,  the  behavior 
of  the  model  during  its  adaptation  is  not  as  well 
understood.  It  is  speculated  that  either  the  process 
will  lock  onto  K  —  2  targets,  ignoring  the  rest,  or 
some  model  components  will  lock  onto  multiple 
targets  sharing  similar  parameters.  In  a  future  study. 


Y,P(WP"n  I  0)  =  &0 (yfjmn)  =  ^Ekm'  ( 10) 

n  n  k 

IV.  PARAMETER  ESTIMATION 

The  goal  is  to  find  the  set  of  model  parameters 
0  =  {Ekm,Qk}  providing  the  best  match  between 
the  the  model  p( wjmn  |  0)  and  the  data  p0( wjmn). 

The  model,  described  in  Section  III,  is  completely 
specified  by  the  the  mixture  weights  Ekm, 
the  variances  ajlk  and  ajk,  and  the  parameters 
{xk,yk,x'k,y'k,xk,yk}  describing  the  target  trajectories. 

In  this  paper  the  Einsteinian  log-likelihood  [22, 
ch.  4.4.1] 

LL(0)  =  ^Po(Wjmn)lnp(w;mn  |  0) 
j,m,n 

K 

=  Y  po(w  jmn  )]nYEk>nPk(Wjmn  I  0Jfc) 
j,m,n  k=  1 

(li) 

will  serve  as  a  criterion  to  quantify  the  similarity 
between  the  model  and  the  data.  Given  the  constraint 
in  (10),  it  can  be  shown  that  LL  will  be  maximized 
when  the  model  and  the  data  match,  i.e.,  for 
piv/jmn  I  0)  =  Po(y?jmn)-  Note  that  maximizing  LL  is 
equivalent  to  minimizing  the  cross-entropy,  a.k.a,  the 
Kullback-Leibler  (KL)  divergence,  a  well-established 
metric  from  the  information  theory  arena  [41]. 
Mathematical  and  physical  justifications  for  both 
Einsteinian  log-likelihood  and  the  KL  divergence 
have  been  discussed  in  detail  [22,  37-40].  It  should 
be  emphasized  that  the  tracking  algorithm  described 
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here  is  a  “batch”  processing  approach,  i.e.,  the  frames 
from  all  time  indices  j  are  processed  jointly.  This  is 
apparent  in  (11)  since  the  log-likelihood  includes  a 
summation  over  j. 

The  mathematical  details  of  the  optimization 
are  contained  in  Appendix  A.  Briefly,  a  system 
of  equations  is  derived  which  is  satisfied  by  the 
parameters  maximizing  LL,  subject  to  the  constraint 
of  (10).  Unfortunately,  an  analytical  solution  is 
intractable  since  the  system  of  equations  is  large, 
coupled,  and  nonlinear.  However,  iterative  techniques 
can  be  employed  to  solve  the  system  of  equations 
and,  in  Appendix  A,  an  efficient  recursive  technique 
is  described  which  is  a  special  case  of  EM.  This 
technique  defines  a  recursive  update  formula  for  each 
parameter.  For  example,  using  the  notation  E^n  and 
0[7)  to  indicate  the  estimates  of  the  parameters  at  the 
7th  iteration,  the  recursive  update  formulas  for  Ekm, 

°dk>  and  are 


pU+i) 
^ km 


(4t)(/+1) 


and 


(4)(/+1) 


where 


jJ2pv(v/j>™)P<r>(k  I  pnn)  (12) 

j,n 

Ej,m,nPo(WjnJP(I)(k  I  jmil)(d jnm  -  D jkJ2 
Po(^jmn)P(l>(k  I  Jmn) 

(13) 

Po(Wjmn)P' [I)(k  I  jmn)(rjmn-Rjk„f 

Y/j,m,nPo(wjmn)P(I)(k\jmn) 

(14) 


P(I\k  |  jmn)  - 


EhnPk^jmn  I  Q^) 
Ek'Ek'lPk'(^jmn  I  Q*?) 

EkmPk^jmn  I  Q^’) 


P^jmn  I  ©(/)) 


(15) 


There  is  an  analogous  rule  for  updating  the  tracking 
parameters  {xk,yk,xk,yk,xk,yk},  which  is  described 
later  in  this  section.  The  equations  above  imply 
starting  with  an  initial  guess  for  the  parameters 
{it®,©®},  then  alternating  between  updating 
P{I)(k  |  jmn)  using  (15),  then  updating  the  parameters 
{£d^1|,0[/+1)}  using  (12)-(14).  In  fact,  these  steps 
correspond  to  the  E-step  and  the  M-step,  respectively, 
of  the  EM  algorithm.  This  iterative  procedure  is 
guaranteed  to  converge  to  a  local  maximum  of  LL,  as 
shown  in  Appendix  B.  Note  that  the  mixture  weights 
Ekm  are  updated  using  (12)  for  all  target  and  clutter 
components  k  =  1,2,..., K.  In  fact,  (12)  is  appropriate 
for  updating  the  weights  of  arbitrary  mixture  models, 
not  only  the  tracking  model  considered  here.  The 
variances  crjk  and  a\ \  are  updated  for  the  appropriate 
model  components,  as  indicated  by  the  target  and 
clutter  model  equations.  Since  the  variance  parameters 


are  initialized  to  large  values,  the  components  are 
initially  fuzzy  and  indistinct.  However,  with  increasing 
iterations,  the  variances  tend  to  adaptively  decrease, 
and  the  individual  components  gradually  converge  and 
lock  on  to  the  target  signatures,  or  to  portions  of  the 
clutter.  Thus,  the  model  iteratively  adapts  to  fit  the 
data,  as  we  demonstrate  with  the  results  in  Section  VI. 

If  the  model  is  viewed  probabilistically 
[22,  ch.  4.4.8],  then  (15)  is  simply  a  form  of 
Bayes’  rule,  and  the  quantities  P,r,(k  \  jmn )  can 
be  construed  as  the  probabilities  that  the  energy 
in  pixel  ( j,m,n )  originates  from  target  or  clutter 
component  k.  Therefore,  Pa>(k  \  jmn )  are  referred 
to  as  the  “association  probabilities”  [22].  Equation 
(12)  therefore  makes  intuitive  sense — it  simply  states 
that  Ek+ 1 ;  is  the  sum  of  all  pixel  values  Pc/w //««)’ 
where  the  pixels  are  weighted  by  their  association 
probabilities.  Similarly,  (13)  and  (14)  correspond 
to  the  usual  definition  of  sample  variance,  with 
the  distinction  that  the  pixels  are  weighted  by  their 
association  probabilities.  From  (15)  it  is  apparent  that 

^2pU)(k  |  jmn)  =  1  (16) 

k 


for  any  iteration  7. 

For  compactness,  we  define  the  angled  bracket 
notation 

(*)(/)  =  J2p0 (y,jmn)pV)(k  I  jmn)(*)  (17) 

j .« 


where  the  asterisk  *  denotes  a  generic  quantity.  Thus, 
( 1 2) — ( 1 4)  can  be  rewritten  in  the  compact  form 


r(/+D  _ 
km 


2  n(/+  1) 


and 


(<*)' 


(l)(/) 

J 

EmiWjnm-Djkm)2)^ 

Em(1)(/) 

y  .  _  f? 

jnm  ikn 


2  \(/+l)  _  E*((rjiB  Pjkm)  ) 


(18) 

(19) 


(20) 


The  converged  values  of  Doppler  variance 
parameters  a^k  are  related  to  the  radar  Doppler 
resolution  and  target  dynamics,  and  are  thus  difficult 
to  estimate  a  priori.  On  the  other  hand,  the  converged 
values  of  the  range  variance  parameters  are 
related  to  the  radar  range  resolution,  modified  by 
the  range  extent  of  the  target,  and  are  assumed  to 
be  better  known  a  priori.  Therefore,  rather  than 
adaptively  estimating  range  variance  using  (20),  it 
may  be  advantageous  to  evolve  ark  according  to  a 
predetermined  schedule  [22,  25,  26].  For  example, 
ark  can  be  initialized  to  a  rather  large  value,  so  that 
the  target  components  can  “see”  large  sections  of  the 
pixel  data,  then  decreased  according  to  an  exponential 
decay  during  EM  iterations  to  a  steady  state  value 
corresponding  to  the  appropriate  sensor  resolution. 
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Finally,  we  derive  the  recursive  update  formula  for 
each  of  the  tracking  parameters  {xk,yk,x'k,y'k,xk,y'k}, 
k  =  1,2,3,. ..,(K  —  2).  In  the  single-sensor  case, 
described  below  in  Section  V,  the  update  formula  for 
these  parameters  is  a  simple  closed-form  expression, 
analogous  to  (18)-(20).  Unfortunately,  a  closed-form 
update  formula  cannot  be  derived  for  the  general 
multi-sensor  case.  However,  local  convergence  is 
still  guaranteed  (see  Appendix  B)  if  the  tracking 
parameters  are  simply  “nudged”  along  the  uphill 
gradient  of  LL  during  each  iteration.  This  actually 
corresponds  to  a  special  case  of  the  generalized 
EM  (GEM)  procedure  [11].  Using  sk  as  a  generic 
placeholder  for  any  one  of  the  tracking  parameters 
{xk,yk,x'k,y'k,x'j.,y'k},  the  update  formula  is  given  by 


„(/+!)  _  M) 


+  h  ■ 


dLL 

dsk 


where  h  is  the  gradient  ascent  stepsize,  and 


(21) 


dLL  _  /  /  rjnm  Pjkm  \  f  jkm  \  \ 

^  yyi  '  \  ® flc  )  V  ^  ^k  )  / 


This  expression  is  derived  in  Appendix  A.  It  is 
straightforward  to  compute  explicit  expressions  for 
dRjkm/dsk  and  ODjkm/dsk  for  each  of  the  parameters 

sk  =  {xk^k^y'k^y'k}  usin§  (3)  and  (4)-  For 
example,  dRjkm/dx°k  =  ~{xm  -xjk)/Rjkm.  Although 
(21)  describes  a  single  gradient  ascent  step  during 
each  iteration,  multiple  steps  can  be  taken  during  each 
iteration  while  holding  the  association  probabilities 
P(k  |  jmn)  fixed.  In  practice,  the  convergence  rate 
can  be  optimized,  using  trial-and-error  to  determine  a 
suitable  combination  of  stepsize  lr  and  steps/iteration. 

Although  the  present  study  is  mainly  concerned 
with  tracking,  it  is  important  to  consider  the  issue 
of  automatic  detection  and  target  declaration.  The 
standard  detection  approach,  e.g.,  described  in  [22, 
sect.  7.2.9],  utilizes  a  log-likelihood  ratio  test  which 
operates  on  the  converged  parameter  values.  The 
issue  of  detection  will  be  studied  in  more  detail  in  the 
future. 


V.  RANGE-ONLY  TRACKING  FROM  A  SINGLE 
SENSOR 


are  still  updated  using  ( 1 8)— (20).  However,  in  the 
single-sensor  case  the  target  track  parameters  can  be 
updated  using  a  closed-form  formula,  which  is  simpler 
and  more  efficient  than  the  gradient  ascent  procedure 
described  by  (21)  used  for  the  multi-sensor  case.  Note 
that  the  analysis  in  this  section  is  very  similar  to  the 
analysis  given  in  [22,  ch.  7.2.6].  However,  whereas  the 
the  previous  work  utilized  a  constant  velocity  model, 
here  a  constant  acceleration  model  is  used  which  is 
slightly  more  complicated. 

Suppose  the  range  between  the  targets  and  the 
(single)  sensor  can  be  described  using  a  constant 
acceleration  model,  i.e., 

Rjkm  =  R°k+R'ktj+R'tf  (23) 


where  Rk  is  the  time-zero  range  of  target  k,  Rk  is  the 
time-zero  range-rate,  and  Rk  is  proportional  to  the 
range-acceleration.  Although  there  is  only  a  single 
sensor  m  =  1,  the  m  subscript  is  maintained  in  the 
quantity  Rjkm  for  consistency  with  the  expression  for 
the  target  model  in  (5).  The  Doppler  (range-rate)  at 
time  tj  is  the  derivative  of  range  with  respect  to  time 
evaluated  at  t •  ,  i.e., 

Djkm=K  +  2Ktj-  (24) 


In  Appendix  A,  an  update  formula  is  derived  for 
iteratively  adjusting  the  values  of  {Rk,Rk,Rk}  in  order 
to  maximize  LL  which,  as  in  the  multi-sensor  case,  is 
described  by  (11).  This  update  formula  is  described 
by  the  following  3  x  3  set  of  linear  equations: 


H, 
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Here,  for  example,  (7?®)(/+1)  denotes  the  updated 
estimate  for  Rk  computed  at  the  (I  +  l)th  iteration.  For 
the  sake  of  compactness,  the  (/)  superscript  has  been 
omitted  from  the  angled  bracket  notation  {*)  used  for 
the  elements  in  the  above  matrices.  However  it  should 
be  understood  that  here 


(*)  =  Po(^jmn)P<l>(k  I  jmn'K *)•  (26) 

j,n 


If  only  a  single  sensor  platform  is  available,  and 
measurements  of  azimuth  are  unavailable,  it  may  still 
be  desirable  to  perform  range-Doppler-only  tracking. 
Here,  the  target  and  interference  (clutter  plus  noise) 
models  given  by  (5) — (7)  remain  valid,  as  well  as 
the  total  model  for  the  data  given  by  (8).  Also,  the 
variance  (a\ \k,cr;.k )  and  mixture  weight  Ekm  parameters 


It  can  be  shown  that,  in  the  absence  of  Doppler 
measurements,  the  set  of  equations  (25)  is  equivalent 
to  second-order  polynomial  regression  in  which 
data  samples  are  weighted  by  their  association 
probabilities. 

Unlike  the  update  formulas  for  Ekm,  <rdk,  and 
ark,  given  by  (18)-(20),  the  matrix  expression 
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Fig.  2.  Upper  left:  scan  frame  data  acquired  at  time  t  =  0  s  in  which  a  single  target  signature  appears  at  Doppler  and  range  values  of 
roughly  (—12  m/s)  and  (258  m),  respectively.  Remaining  three  plots  show  evolution  of  model,  which  adapts  to  fit  data  as  iterations 
increase.  White  Xs  in  lower  right  plot  indicate  converged  range/Doppler  coordinates  R®  and  R'k  for  target  components  in  model. 


(25)  actually  represents  a  set  of  three  coupled 
update  formulas  for  the  three  tracking  parameters 
{Rk,R'k,R'l}.  Upon  each  iteration,  these  parameters 
are  updated  by  inverting  (25).  Note  that  a  separate 
3x3  matrix  inversion  is  performed  for  each  target 
component  A  =  1,2,. ..,(K  —  2). 

VI.  RESULTS 

In  this  section,  we  describe  results  when  the 
algorithm  is  applied  to  two  cases.  First,  the  algorithm 
is  tested  against  experimental  stretch  radar  data 
[33]  from  a  single  sensor  using  only  the  range  and 
Doppler  data  from  the  sensor.  The  single  sensor 
results  are  presented  as  a  proof  of  concept  and  to 
illustrate  the  performance  of  the  algorithm  in  realistic, 
inhomogeneous  clutter,  variable  signal  to  interference, 
and  accelerating  and  maneuvering  targets.  In  the 
second  case,  we  consider  multiple  sensors  viewing 
the  same  scene.  The  sensors  are  assumed  to  have  poor 
or  no  azimuth  resolution  but  good  range  and  Doppler 
resolution.  Combining  the  data  from  multiple  sensors 
offers  the  opportunity  to  more  accurately  locate 
and  track  targets  using  triangulation  of  the  sensor 
data  if  the  difficult  problem  of  target  association 
can  be  solved.  In  the  second  case  we  present  results 
from  computational  experiments  designed  to  test  the 
algorithm  since  experimental  multi-sensor  data  are  not 
available  to  us. 

Let  us  now  describe  the  single-sensor  results 
computed  from  experimental  data,  using  the  simplified 
version  of  the  algorithm  described  in  Section  V.  These 
results  will  be  mainly  qualitative  since  neither  true 


target  tracks  nor  detailed  information  about  the  radar 
were  provided.  Nevertheless,  this  analysis  serves 
as  a  useful  proof-of-concept  for  the  algorithm.  The 
experimental  setup  consisted  of  a  stretch  receiver 
mounted  on  a  tower  overlooking  a  wooded  area, 
and  data  were  collected  as  various  targets  moved 
through  the  area.  The  raw  data  were  converted  to  a 
suite  of  range/Doppler  scan  frames  using  standard 
processing  techniques  [33],  as  discussed  in  Section  II. 
The  algorithm  described  in  Section  V  was  used  to 
jointly  process  multiple  sets  of  scan  frames,  where 
each  scan  frame  is  computed  from  data  centered 
around  a  different  time  f  •,  j  =  1,2,..., J. 

The  upper  left  plot  of  Fig.  2  shows  a  single  scan 
frame  acquired  at  a  reference  time  of  /  =  0  s.  A 
target  signature  appears  as  a  “blob  of  energy”  at 
Doppler  and  range  values  of  roughly  (—12  m/s)  and 
(258  m),  respectively.  There  is  a  significant  ridge 
of  ground  clutter  from  the  stationary  background 
centered  at  zero  Doppler,  extending  across  all  range 
bins.  There  is  also  a  noise  background  which  is 
roughly  uniform  in  range/Doppler,  although  with 
some  significant  inhomogeneities.  For  the  model 
we  used  five  components:  three  target  components, 
one  uniform  background  noise  component,  plus  one 
clutter  component  shaped  like  a  ridge  centered  at  zero 
Doppler. 

The  remaining  three  plots  of  Fig.  2  show  how 
the  model  evolves  with  increasing  iterations.  Since 
the  model  variances  are  initialized  to  large  values, 
the  model  is  initially  broad,  fuzzy,  and  indistinct 
(upper  right).  However,  as  the  iterations  progress,  the 
parameters  adapt  in  such  a  way  as  to  eventually  define 
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Fig.  3.  Same  as  Fig.  2,  however  here  scan  frame  was  acquired  at 

and  range  values  of  roughly  (— ! 
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a  later  time  of  /  =  6  s.  Target  signature  has  now  moved  to  Doppler 
i  m/s)  and  (195  m),  respectively. 


Data,  time  =  6  seconds  Model,  iteration=  300 
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Fig.  4.  Left-hand  plot  shows  same  scan  frame  shown  in  Fig.  3,  however  a  large  amount  of  artificial  noise  was  added  which  almost 
completely  obscures  target  signature.  Nevertheless  (right-hand  plot),  one  of  the  model  components  was  able  to  adapt  to  range/Doppler 

coordinate  in  close  vicinity  of  target  signature. 


a  model  which  shares  most  of  the  significant  features 
of  the  data.  The  lower  right  plot  in  Fig.  2  shows 
the  converged  model  after  200  iterations.  The  white 
Xs  indicate  the  converged  range/Doppler  estimates 
R J?  and  R'k  for  the  three  target  components.  Notice 
that  one  of  the  three  target  components  has  properly 
converged  to  the  true  target  signature,  while  the  extra 
two  target  components  have  converged  to  portions 
of  the  inhomogeneous  background  clutter.  Although 
Fig.  2  shows  the  data  and  model  corresponding  to 
only  a  single  scan  frame,  it  is  important  to  note  that 
the  results  were  computed  by  jointly  processing 
10  scan  frames,  starting  at  t  =  0  and  proceeding  in 
intervals  of  0.25  s. 

The  upper  left  plot  of  Fig.  3  shows  scan  frame 
data  corresponding  to  a  different  time,  t  =  6  s.  The 
target  signature  has  now  moved  to  Doppler  and  range 
values  of  roughly  (—9  m/s)  and  (195  m),  respectively. 
Here,  the  inhomogeneities  in  the  background  are 
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even  more  extreme  than  for  the  t  =  0  case  shown  in 
Fig.  2.  However,  the  model  was  able  to  adapt  quite 
nicely,  with  one  of  the  target  components  locking  on 
to  the  true  target  signature,  and  the  extra  components 
adapting  to  fit  the  clutter. 

In  Figs.  2  and  3,  the  target  signatures  were  rather 
strong,  with  signal-to-interference  ratios  of  roughly 
3  dB  (computed  from  the  peak  value  of  the  signature 
relative  to  an  average  background  value  in  the  image 
frame).  Next,  the  tracking  experiment  was  repeated 
for  the  same  case  shown  in  Fig.  3.  However  the 
problem  was  intentionally  made  more  difficult  by 
adding  a  large  amount  of  artificial  noise,  so  that  the 
signal-to-interference  ratio  was  reduced  to  roughly 
0  dB.  This  noise  is  uniformly  distributed  in  range 
and  Doppler,  and  is  signal-independent,  i.e.,  it  was 
simply  added  to  the  intensity  image.  The  data  and 
the  converged  model  for  this  noisy  experiment  are 
shown  in  Fig.  4,  which  shows  that  one  of  the  model 
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Fig.  5.  Tracking  over  extended  time  interval  for  smoothly  decelerating  target.  Each  scan  frame  corresponds  to  different  time  instance 
with  12  s  interval.  Overlaid  is  dark  line  showing  track  computed  using  our  method. 


components  locked  on  to  the  target  signature  in  the 
data,  despite  the  fact  that  the  target  signature  is  almost 
completely  obscured  by  noise. 

Although  the  present  study  is  mainly  concerned 
with  tracking,  it  is  important  to  consider  the  issue 
of  automatic  detection  and  target  declaration.  For 
example,  in  the  converged  model  shown  in  the  lower 
right  plot  of  Fig.  3  there  are  three  target  components 
which  converged  to  the  three  different  range/Doppler 
coordinates  indicated  with  white  Xs.  One  component 
has  locked  on  to  the  target  signature,  while  the  other 
two  have  adapted  to  fit  portions  of  the  clutter.  Here, 
the  detection  problem  consists  of  automatically 
declaring  the  true  target,  while  disregarding  the  other 
two  components. 

The  standard  detection  approach,  e.g.,  described 
in  reference  [22,  sect.  7.2.9],  utilizes  a  log-likelihood 
ratio  test  evaluated  at  the  converged  track  and  variance 
parameter  values.  In  the  examples  shown  thus  far,  it 
appears  that  perhaps  a  simple  detection  rule  could 
be  devised  based  upon  the  compactness  of  the  target 
component,  i.e.,  based  upon  the  converged  values 
of  the  variance  parameters  adk  and  a rk.  The  issue  of 
detection  will  be  studied  in  more  detail  in  a  future 
publication. 

Next,  sample  results  are  presented  in  which 
tracking  was  performed  over  longer  time  intervals. 
Here,  a  sliding  window  approach  was  used  to 
estimate  tracks  along  extended,  irregular,  paths  in  the 
range/Doppler  space.  Overlapping  sets  (“batches”)  of 
6  scan  frames  were  used  to  compute  each  point  along 
the  path.  For  example,  the  first  point  on  the  path  was 
computed  by  jointly  processing  6  scan  frames  within 
the  time  interval  0  <tj<\  .25  s,  the  second  point 


on  the  path  was  computed  by  jointly  processing  6 
scan  frames  within  the  time  interval  0.5  <  t  ■  <  1.75  s, 
etc.  Fig.  5  shows  results  from  tracking  a  smoothly 
decelerating  target.  Four  different  scan  frames  are 
shown  corresponding  to  times  of  t  =  0,  2.5,  6,  and 
12  s.  Superimposed  on  the  images  is  a  dark  line 
which  indicates  the  extended  track  computed  using 
the  sliding  window  method  over  a  12  s  interval. 
Notice  the  close  agreement  between  the  location  of 
the  target  signatures  and  the  computed  track.  Next, 
Fig.  6  shows  results  from  tracking  a  maneuvering 
target.  Here  the  target  accelerates  significantly  over 
the  time  observation  interval,  whereas  previously 
the  target  velocity  changed  relatively  slowly.  Each 
scan  frame  in  Fig.  6  corresponds  to  a  different  time 
instance  within  a  15.5  s  interval.  Overlaid  is  a  dark 
line  showing  the  track  computed  using  the  sliding 
window  method  which,  again,  shows  close  agreement 
between  the  location  of  the  target  signatures  and  the 
computed  track. 

In  the  remainder  of  this  section  results  are 
presented  for  the  full,  multi-sensor,  version  of  the 
algorithm  described  in  Section  IV.  Since  experimental 
multi-monostatic  data  were  unavailable  to  us,  the 
results  are  based  on  sets  of  synthetic  range/Doppler 
scan  frames  of  the  type  that  would  be  produced  from 
UHF  stretch  radar  data  [33].  In  these  simulations, 
emphasis  was  more  on  data  association  than  on 
detection  performance.  Therefore,  the  data  were 
generated  with  a  relatively  high  signal-to-interference 
ratio  of  approximately  8  dB,  computed  from  the 
peak  value  of  the  signature  relative  to  an  average 
background  value  in  the  image  frame.  However, 
the  zero-Doppler  ridge  from  stationary  clutter  was 
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Fig.  6.  Tracking  over  extended  time  interval  for  maneuvering  target.  Each  scan  frame  corresponds  to  different  time  instance  within 
15.5  s  interval.  Overlaid  is  dark  line  showing  track  computed  using  our  method. 
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Fig.  7.  Geometry  for  multi-sensor  simulation.  Each  arrow  shows 
direction  of  motion  for  corresponding  target,  and  numbers  in 
parentheses  are  x  (east)  and  v  (north)  components  of  velocity 

(44)- 

generated  with  roughly  the  same  peak  amplitude  as 
the  peak  amplitude  of  the  target  signatures.  Therefore 
this  ridge  overlaps  and  obscures  the  target  signatures 
having  low  Doppler,  as  illustrated  in  Fig.  8.  The 
synthetic  data  were  generated  by  placing  “blobs”  of 
energy  at  appropriate  target  and  clutter  locations  at 
each  frame  time  without  errors  or  biases,  other  that 
those  associated  with  the  dimensions  of  the  “blobs” 
which  are  related  to  system  resolution  and  target 
spreading  effects.  While  errors  and  biases  may  be 
important  in  the  ultimate  application,  we  want  to 
illustrate  by  this  example  the  native  data  association 
performance  of  the  proposed  multi-sensor  algorithm 
with  multiple  targets  and  sensors. 

Fig.  7  shows  the  geometry  for  the  experiment. 
There  are  5  sensors  and  4  targets,  two  of  which 


have  constant  velocity  and  two  of  which  are  turning 
and  accelerating.  These  four  targets  have  x  (east) 
and  y  (north)  components  of  their  velocities  given 
by  (x'k,y'k)  =  (0,-4), (7, 5), (3, 2), (-15,1)  m/s,  and 
corresponding  acceleration  components  proportional 
to  (4',y")  =  (-0.4,1), (-2,0), (0,0), (0,0)  m/s2.  Sensors 
1-3  are  spaced  roughly  10  meters  apart,2  while  the 
positions  of  sensors  4  and  5  are  more  widely  spaced. 
The  spatial  diversity  is  exploited  to  “triangulate” 
accurate  track  estimates  for  the  targets.  The  targets 
were  observed  by  all  five  sensors  at  4  different  times 
separated  by  0.25  s.  During  the  4  frames  of  data  the 
targets  move  a  distance  comparable  to  the  sensor 
range  resolution.  Thus,  the  algorithm  was  required 
to  jointly  process  a  total  of  twenty  scan  frames  in 
order  to  estimate  target  tracks.  Fig.  8  shows  four  of 
these  twenty  scan  frames,  corresponding  to  sensors 

1  and  5  at  times  tx  =  0  and  t4  =  0.75  s.  From  this 
figure  it  is  apparent  that  the  data  association  problem 
is  not  trivial.  For  example,  it  is  not  obvious  “by  eye” 
which  signatures  in  the  bottom-right  plot  of  Fig.  8 
correspond  to  which  signatures  in  the  top-left  plot. 

The  model  was  constructed  of  8  components,  i.e., 

6  target  components  plus  a  uniform  noise  component, 
plus  a  zero-Doppler  clutter  ridge  component.  We 
assumed  no  prior  knowledge  about  how  many  targets 
were  actually  present,  and  therefore  2  more  target 
components  were  used  than  actual  targets  present. 

The  code  seemed  to  converge  most  quickly  by  starting 
off  with  30  iterations  using  data  only  from  sensors 
1-3  (these  sensors  are  relatively  close  together), 

2  Sensors  1-3  do  not  operate  as  an  array  since  their  data  are 
processed  noncoherently. 
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Fig.  8.  Synthetic  range/Doppler  image  data  showing  scan  frames  for  sensor  1  (time  =  0,0.75  s)  and  sensor  5  (time  =  0,0.75  s).  There 
are  4  target  signatures  present.  Note  that  signatures  of  lower  velocity  targets  significantly  overlap  clutter  ridge  centered  at  zero  Doppler. 
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Fig.  9.  Model  evolution  for  sensor  1.  time  =  0  data.  Compare  converged  model  (bottom-right  plot  here)  to  actual  data 

(top-left  plot  of  Fig.  8). 

followed  by  80  additional  iterations  after  adding  the  of  the  model.  Separate  plots  are  given  for  the  east 
data  from  sensors  4  and  5.  Note  that  for  each  of  these  (x)  components  of  time-zero  velocity  x'k,  time-zero 

iterations,  the  tracking  parameters  were  adjusted  using  position  xk,  and  acceleration  x!'k.  Similar  plots  were 

50  gradient  ascent  steps,  as  discussed  in  Section  V.  obtained  for  the  y  (north)  components  y'k,  yk .  and  y". 

The  evolution  of  the  model  is  shown  in  Fig.  9  but  they  are  not  shown  here.  In  Fig.  10,  a  plot  is  also 

for  the  scan  frame  acquired  by  sensor  1  at  t]  -0  s.  given  showing  the  evolution  of  the  mixture  weights 

The  model  is  initially  fuzzy  and  indistinct  (top-left  Ek  for  each  of  the  target  components  (Ek  of  clutter 

plot),  but  after  30  or  more  iterations  it  converges  to  an  components  are  not  shown). 

accurate  approximation  of  the  data,  as  can  be  seen  As  mentioned  above,  the  first  30  iterations  were 

by  comparing  the  bottom-right  plot  in  Fig.  9  with  performed  using  data  only  from  sensors  1-3  while 

the  top-left  plot  in  Fig.  8.  In  Fig.  10  it  is  shown  how  the  remaining  iterations  were  performed  using  all 

the  parameter  values  evolve  during  the  adaptation  data  from  sensors  1-5.  There  are  curves  in  the 
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xk’  -  time-zero  velocity,  east  (m/s)  x°  =  time-zero  position,  east  (m) 


xk”  =  acceleration,  east  (m  1^) 
1  r 


iteration 


mixture  weight  Ek 


iteration 


Fig.  10.  Evolution  of  parameter  estimates  with  iterations.  Notations  cl,  c2,  etc.,  in  legend  refer  to  target  components  1-6  used  in 
model.  Although  we  only  show  plots  for  x  (east)  components  of  position,  velocity,  and  acceleration,  similar  plots  were  obtained 

for  y  (north)  components. 


plots  of  Fig.  10  for  each  of  the  6  target  components 
labelled  “cl-c6,”  and  these  are  compared  with 
the  true  parameter  values  (dashed  lines).  Recall, 
there  are  only  4  targets  actually  present.  From  the 
plots,  it  can  be  seen  that  components  cl,  c4,  and  c5, 
converge  asymptotically  to  the  true  target  parameters. 
In  contrast,  component  c2  seems  to  drift  around 
without  locking  onto  a  target,  and  its  corresponding 
mixture  weight  parameter  decreases  toward  zero 
with  increasing  iterations.  Based  upon  its  dwindling 
Ek  value,  component  c2  could  easily  be  pruned 
from  the  model.  The  most  interesting  case  involves 
components  c3  and  c6,  which  seem  to  lock  onto  the 
same  target,  as  shown  by  the  fact  that  their  velocity, 
position,  and  acceleration  estimates  converge  together 
in  order  to  share  a  true  target  track.  The  mixture 
weight  parameter  Ek  plot  (right-most  plot)  shows 
that  the  weights  for  c3  and  c6  roughly  sum  up  to  the 
true  mixture  weight  of  100.  Presumably,  it  would  be 
straightforward  to  detect  this  condition  in  which  two 
components  lock  onto  the  same  target.  Then  one  of 
the  two  could  simply  be  pruned  from  the  model  if 
necessary.  It  is  an  open  question  how  the  model  would 
adapt  for  cases  in  which  the  true  number  of  targets 
exceeds  the  number  of  model  components.  Strategies 
for  handling  this  so-called  “model  selection”  problem 
for  general  mixture  models  are  described  in  [34],  [35], 
[36]  and  the  references  therein. 

VII.  CONCLUSIONS 

In  this  paper  we  describe  a  mathematical  algorithm 
for  robust  multi-target  tracking  from  multiple  radar 


platforms,  for  difficult  cases  in  which  measurements 
of  azimuth  are  unavailable.  A  simplified  version  of  the 
algorithm  is  also  presented  which  is  appropriate  for 
tracking  from  a  single  radar  platform. 

The  single-sensor  version  of  the  algorithm  is  tested 
on  experimental  data.  These  results  show  the  approach 
to  be  very  promising,  and  robust  in  the  presence  of 
significant,  inhomogeneous,  background  clutter.  The 
full,  multi-sensor,  algorithm  is  tested  on  synthetic 
data.  These  results  demonstrate  that  accurate  tracks 
can  be  estimated  by  exploiting  spatial  diversity  in  the 
sensor  locations.  Furthermore,  the  algorithm  appears 
to  be  robust  in  the  presence  of  significant  clutter,  and 
uncertain  knowledge  regarding  the  number  of  targets 
present. 

It  should  be  emphasized  that  while  the  initial 
results  are  promising,  they  are  only  preliminary. 
Significant  issues  require  further  study  in  order  to 
make  the  algorithm  practical.  For  example,  it  is  well 
known  that  EM  can  converge  to  local  maxima  in  the 
optimization  criterion,  and  therefore  initialization 
will  likely  affect  the  converged  solution.  We  still 
need  to  determine  how  often  these  “traps”  will  occur, 
under  which  conditions  they  can  cause  the  algorithm 
to  fail,  and  how  we  can  initialize  the  algorithm  to 
prevent  them.  Also,  since  the  true  number  of  targets 
is  seldom  known  a  priori,  more  work  is  needed  to 
develop  strategies  for  choosing  the  number  of  target 
components,  pruning  unneeded  components,  and 
adding  components  when  necessary.  Furthermore, 
automatic  detection  strategies  and  performance  will 
need  to  be  studied  to  judge  our  ability  to  distinguish 
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targets  from  nonhomogeneous  clutter.  Finally,  more 
study  is  required  to  determine  the  robustness  of  the 
algorithm  in  the  presence  of  sensor  platform  errors 
and  vibrations. 


seek  to  minimize  the  Lagrangian 

F  =  -LL  +  A  (  -  J2Po(wjmn)  )  (30) 

\  k'  n  ) 


APPENDIX  A.  DERIVATION  OF  PARAMETER 
ESTIMATION  EQUATIONS 

Here  the  equations  are  derived  which  are  satisfied 
by  the  parameters  maximizing  LL  subject  to  a 
constraint.  The  general  procedure  is  first  developed, 
valid  for  an  arbitrary  mixture  model.  These  results 
are  then  used  to  derive  the  particular  equations 
corresponding  to  the  tracking  model.  Good  references 
for  the  general  procedure  described  here  are 
Perlovsky’s  book  [22],  and  the  seminal  work  by  Duda 
and  Hart  [32].  The  notation  we  use  is  very  similar  to 
[32], 

First,  consider  maximizing  LL  using  a  general 
mixture  model.  One  might  seek  to  maximize  LL 
directly,  by  finding  values  for  the  set  of  parameters 
Qk  for  which  all  partial  derivatives  8LL/dOk  are  zero. 
Using  the  expression  for  LL  from  (11), 


with  respect  to  Ekm,  where  A  is  the  Lagrange 
multiplier.  Setting  the  partial  derivative  of  F  to  zero, 

dF  dLL 

ofT  ~~dFT  +  x-°  (31) 

UCjkm  UCjkm 

the  following  value  for  A  is  obtained 

dLL 


A  = 


dE, 


(32) 


km 


An  expression  for  dLL/8Ekm  is  now  needed,  and  its 
derivation  is  similar  to  (27).  We  obtain 

1  J^g^-^kmPk^jmn  I  0*)1 

km  j,n  km 


Combining  this  with  (32), 


dLL 

dek 


j,m,n 


]n'y2Ek'r„Pk'(w  jmn  I©*-) 
k' 


J  N 

XEkm  = 

j  =  1  «=1 


(33) 


=  £ 

j,m,n 


Po^jmn) 

_^2k'  Ek'mPk'(yl jmn  I  ©*'). 


Q 0  I  E k/nP I  ©A:  1 1 


=  'yZPo^jmn) 

j,m,n 


EkmPk(W, jmn  I  ©*) 

_^2k'EklmPk'('IVjmn  I  ©t')_ 


Next  we  perform  the  summation  over  k  on  both  sides 
of  (33),  then  use  (10)  and  (16)  to  simplify.  This 
results  in  the  solution  for  the  Lagrange  multiplier 
A  =  J,  where  J  is  the  number  of  scan  frames  (time 
indices)  from  each  sensor  used  in  the  estimation. 
Substituting  this  value  back  into  (33),  we  obtain 


d 

QQ-\Ek,nPk('»jnn  I  ©*)] 

X  k _ 

EkmPk^jmn  I  ©*) 

=  '52Po('"jmn)P(k  I  imn)-^-MEkmPh(ykJmn  |  0Q1. 

j,m,n  k 

(27) 

Here,  we  made  use  of  the  well-known  identity 
d \ny(x)/dx  =  [  I  /y(x)\c)y(x)/dx,  and  the  definition  of 


P(k  |  jmn)  = 


EhnPk^jnm  I  ©^) 

^2/k'  Ek'mPk'^y^ jmn  I  ®k') 


(28) 


Using  (27),  the  set  of  parameters  Qk  that  maximize  LL 
will  satisfy 

qq-  =  J2Po(wjmn)P(k  I  jmn)gQ-lnPk(wjmn  I  ®k)  = 

k  j,m,n  k 


k  =  1,2 


(29) 


The  maximization  of  LL  must  also  be  performed 
with  respect  to  the  mixture  weights  Ekm.  However  it  is 
now  necessary  to  incorporate  the  constraint  from  (10) 
using  the  method  of  Lagrange  multipliers.  Thus  we 


Ekm  =  jYl  P0(wjmn)p(k  I  jmn)-  (34) 

j,n 

Equation  (34)  is  satisfied  by  the  set  of  mixture 
weights  Ekm  that  maximize  LL ,  subject  to  the 
constraint  in  (10). 

A  set  of  equations  has  now  been  derived  governing 
the  complete  set  of  model  parameters  0  =  {Ekm,Qk}, 
for  k  =  1,2 and  m  =  1,2 Equation  (29) 
governs  the  parameters  Qk,  while  (34)  governs  the 
mixture  weights  Ekm.  Unfortunately,  solving  for 
the  parameters  directly  from  these  equations  can  be 
problematic.  The  difficulty  arises  because  P(k  \  jmn )  is 
a  function  of  all  of  the  unknown  parameters  {Ekm,Qk} 
(see  (28)).  Thus,  if  there  are  Q  unknown  parameters, 
it  would  be  required  to  solve  a  coupled  set  of  Q 
nonlinear  equations.  However,  the  problem  is  not 
hopeless,  and  the  maximum  of  LL  can  be  reached 
in  an  iterative  fashion,  alternating  between  updating 
P(k  |  jmn),  then  updating  the  set  of  parameters 
{Ekm,Qk}.  We  use  the  notation  PaHk  \  jmn )  and 
{Ekm’®k}  t0  indicate  the  values  of  these  quantities 
at  the  /th  iteration.  Then,  the  explicit  update  formula 
for  pV\k  |  jmn )  is  given  by  (15),  while  the  parameters 
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are  updated  in  order  to  satisfy  modified  versions  of 
(29)  and  (34),  i.e., 


and 


Ekm  ' '  =  jJ2P0^jmn)pU)(k  I  Jmn) 

j,n 


(35) 


^2P0^jmn)pV\k  I  imn ) 

j,m,n 


A'nPk(^ 


I©*) 


=  0. 


(36) 

This  iterative  procedure  is  guaranteed  to  converge  to 
a  local  maximum  of  LL,  as  shown  in  Appendix  B.  It 
will  be  of  later  use  to  note  that 


=  T'Poty'jmn)  W 

k  n 


for  all  iterations  I.  This  follows  from  (35)  by  using 
the  fact  that  jTk  P{h(k  \  jmn )  =  1,  as  can  be  shown 
from  (15). 

The  update  formula  in  (35)  for  the  mixture  weight 
parameter  is  a  simple,  closed-form  expression.  The 
precise  form  of  the  update  equation  for  Qk  depends 
upon  the  particular  model  pk(y/:mn  \  Qk)  that  gets 
plugged  into  (36).  In  some  cases,  this  leads  to 
closed-form  expressions,  for  example  when  estimating 
the  covariance  and  mean  parameters  of  Gaussian 
mixtures.  Also,  in  this  paper  closed-form  equations 
are  obtained  when  estimating  adk  using  (13)  or  when 
estimating  the  single-sensor  range-only  tracking 
parameters  using  (25).  However,  sometimes  the  model 
is  specified  such  that  a  closed-form  expression  is 
not  possible,  for  example  when  estimating  the  track 
parameters  using  multi-sensor  data,  as  described  at  the 
end  of  Section  IV.  In  this  case,  as  an  alternative  to  the 
closed-form  update  formula  of  (36),  the  parameters  (~)k 
can  be  updated  by  moving  along  the  uphill  gradient  of 
LL,  i.e., 


0(/+|)  =  Q({>  +  li- 


OH 


<90, 


e*=e® 


=  Gk’  +  h  ■  E  Po(wjmn)pV\k  I  jmn) 

j,m,n 


QQ-klnPk(Wjmn  I  0i) 


©t=e!'> 


(38) 


Here  h  is  the  stepsize  which  can  be  chosen,  for 
example,  by  trial  and  error.  Note  that  the  procedure 
described  by  (38)  is  not  pure  gradient  ascent,  since 
this  step  is  alternated  with  updating  PU)(k  \  jmn )  via 
(15)  and  Ek^  using  (35).  In  fact,  this  procedure  is  a 
special  case  of  GEM  [11].  Note  that  in  some  cases 
the  set  of  parameters  <dk  may  be  split,  with  some 
members  being  updated  using  gradient  ascent  (38) 
and  some  being  updated  using  a  closed-form  equation 
(36). 


The  general  procedure,  developed  above,  is  valid 
for  an  arbitrary  mixture  model.  These  results  are 
now  used  to  derive  the  specific  parameter  estimation 
equations  for  the  tracking  application  considered  in 
this  paper. 

Variance  Parameters :  The  recursive  update 
formula  is  now  derived  for  the  variance  parameter  adk 
In  order  to  use  the  update  formula  given  by  (36),  it 
is  necessary  to  compute  an  explicit  formula  for  the 
factor  within  the  square  brackets,  by  substituting  the 
specific  model  pjfc(w  mn  |  0, )  used  for  the  tracking 
application.  Both  the  target  model  of  (5)  and 
the  clutter  model  of  (7)  give  the  following  result 
(neglecting  irrelevant  terms): 


d 


dGl°dk) 


ln/L(w/m„  I  ©,) 


dG/°d0  l 


yMlA4)-y(l  /^djnm-DjkJ 


2^ (Ik  2  jnm  Ejkm  ) ' 


Note,  it  is  actually  more  convenient  to  work  with  the 
inverse  variance  1  /ajlk  rather  than  the  variance  itself, 
as  shown  in  the  above  equation.  Substituting  this 
expression  into  (36)  leads  to  the  update  formula  for 
the  variance  parameter  described  by  (13).  The  update 
formula  in  (14)  for  the  range  variance  a2k  is  derived  in 
a  similar  fashion. 

Track  Parameters,  Multi-Sensor  Case:  Due 
to  the  model  complexity,  a  closed-form  update 
formula  cannot  be  derived  for  the  tracking  parameters 
in  the  multi-sensor  scenario.  However,  a  GEM 
update  formula  can  be  derived,  starting  with  (38), 
and  computing  the  factor  in  square  brackets  by 
inserting  the  track  model  given  by  (5).  Using  sk  as  a 
generic  placeholder  for  any  of  the  track  parameters 
{xk,xk,xk,yk,yk,yk},  we  compute  (neglecting  irrelevant 
terms) 


d 

|  0,) 


- (V 

9  2  v  jnm 
A(Jrk 


R jkm  ) 


0^-2  (djnm 
A(Jdk 
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jkm 


(7 
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rk 


9Rjkm\ 

dsk  ) 
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D 


jkm 


<T 


2 

dk 


Inserting  this  expression  into  (38),  the  parameter 
update  formula  described  by  (21)  and  (22)  is  obtained. 

Track  Parameters,  Single  Sensor  Case:  The 
problem  of  range-only  tracking  using  a  single  sensor 
is  discussed  in  Section  V.  The  goal  is  to  find  values 
for  the  track  parameters  {Rk,Rk,Rk},  for  each  target 
component  k  =  1,2,3,..., (K  —  2),  which  maximize 
LL.  In  this  case  the  target  model  given  by  (5)  is  still 
valid.  However  the  equations  relating  {Rk,Rk,Rk} 
to  the  quantities  Rjkm  and  D  -km  are  now  given  by 
(23)  and  (24),  respectively.  Although  there  is  only  a 
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single  sensor  m  =  1,  the  m  subscript  is  maintained  in 
the  quantities  R;km  and  Djkm  for  consistency  with  the 
expression  for  the  target  model  in  (5). 

The  update  formulas  for  {Rk,R'k,Rk},  are  derived 
by  starting  with  the  general  formula  of  (36),  and 
substituting  an  explicit  formula  for  the  factor  in 
square  brackets,  which  is  computed  using  the  target 
model  pk{ wjmn  \  Qk)  in  (5).  Using  sk  as  a  generic 
placeholder  for  any  of  the  parameters  {Rk,Rk,Rk},  the 
square-bracketed  factor  in  (36)  is 


B  ^  =  — 

jmn  l  Ax  Qg 


k  L 


\(\/<J2rk){r 


jnm 


'  Rjkm) 


^(l/Vdk)(djnm-Djkmf 


r  —  R 

'jnm  A'jkm 


(OR 


jkm 


urk 

^ jnm  D jkm 


dk 


\  dsk 


9Djkm 

dsu 


Substituting  this  expression  into  (36),  the  parameter 
update  formula 


(39) 


is  obtained  where,  for  compactness,  we  have  used 
the  bracket  notation  (*)a>  defined  in  (17).  Note  that 
(23)  and  (24)  are  used  to  express  Rjkm  and  D-km  and 
their  derivatives  with  respect  to  sk  =  {Rk,Rk,Rk}  so, 
for  example,  dRjkm/dRQk  =  1  and  dRjkm/dR'k  =  tj.  For 
each  target  component  k,  (39)  is  used  to  produce 
three  different  equations  for  the  three  track  parameters 
sk  =  {Rk,Rk,Rk}.  First,  substituting  sk  =  Rk  into  (39), 
we  obtain 


(R0kf+1\l)(n 


+(R,kf+i)(tjr+(KYM\tjr=(r. 


\u) 


(40) 

Here,  for  example,  (RkYI+l)  denotes  the  updated 
estimate  for  Rk  computed  at  the  (7  +  l)th  iteration. 
Note  that  (40)  is  linear  with  respect  to  the  three 
parameters  Rk,  R'k,  and  Rk.  By  substituting  sk  =  R'k, 
then  sk  =  R",  into  (39)  we  obtain  two  additional 
equations  which  are  also  linear  with  respect  to  Rk, 

R'k,  and  R”.  Thus  there  are  three  equations  that  are 
linear  with  respect  to  three  unknowns  which,  together, 
can  be  expressed  using  matrix  notation  as  shown 
in  (25)  given  in  Section  V.  Note  there  is  a  separate 
3  x  3  set  of  equations  for  each  target  component 
k  =  1 , 2, ....  (77  —  2).  In  each  iteration  the  tracking 
parameters  are  updated  by  inverting  these  matrix 
equations. 


APPENDIX  B.  CONVERGENCE  PROOF 


This  convergence  proof  generally  follows  the  one 
given  in  [22,  ch.  5.6.3],  although  here  we  provide 
some  details  that  were  left  out  in  the  reference.  The 
proof  is  valid  for  the  general  mixture  model,  and  is 
not  restricted  to  the  tracking  model  used  in  this  paper. 

In  order  to  show  convergence  to  a  local  maximum 
of  LL,  we  show  that  LL  never  decreases  from  one 
iteration  to  the  next,  i.e.,  LL((-)U+ 1 ')  —  LL((-)a>)  >  0, 
where  LL(0(/)  is  the  log-likelihood  computed  in  the 
7th  iteration.  From  (11), 

ZZ(0(/+1))  -  LL(6(/)) 

=  X/°(w;™)[ln^(w/™ 1  ©</+1))-lnT(w;™  |  e(/))]. 

j,m,n 

(41) 

It  is  easy  to  see  from  (15)  that  J2kP{I\k  \  jmn )  -  1, 
and  therefore  we  can  insert  this  unity  factor  into  the 
above  expression,  i.e., 

LL(0</+1))-LL(0(/)) 

=  Y,Po(^j,m,)Y,p(I^k  I  ’mn) 

j,m,n  k 

X  nnp(w/m„  I  0,/+i))  —  lnp(wjmn  I  0(/))]. 

(42) 

We  can  further  expand  the  expression  by  noting 
lnp(w/m„  |  0(/))  =  In  (E^pk(Yrjm„  \  ©f))  -  In  P{I\k  \  jmn) 

as  can  be  seen  from  the  definition  of  PU)(k  \  jmn )  in 
(15).  Using  this  fact,  and  rearranging  terms,  we  can 
obtain 


LL(0</+1))  —  LL(0(/)) 


=  jj2p0(wj,nn)l>2pt,>(k  I  7ww)ln 

j,m,n  k 


PfI\k  |  jmn) 
PV+l)(k  |  jmn) 


+  E  Po^jmn)PV\k  |  jmn)ME«:"pk(yvjmn  \  0</+1>)) 

j,k,m,n 


-  E  Po^jmJP^tt  I  j'n'P^KmPk^jmn  I  ©f))’ 

j,k,m,n 

(43) 

Since  J2kPin(k  \  jmn)  =  1  and  PU)(k  \  jnm)  >  0 
for  all  iterations  /,  the  log-sum  inequality  [41]  or 
Jensen’s  inequality  can  be  used  to  prove  the  first  term 
on  the  right-hand  side  of  (43)  is  >  0.  Thus,  in  order 
to  show  convergence  we  need  to  show  the  sum  of 
the  second  and  third  terms  is  >  0,  which  can  be  done 
by  analyzing  the  rules  governing  the  updates  of  the 
parameters  {EhirQk}.  The  rule  governing  the  update 
of  the  weights  Ekm  is  given  by  (35),  which  can  be 
rewritten  as 

Epo(w;m«)p</)(fe  I  jmn)— -7=0.  (44) 

j,n  ^ km 
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Then,  since 

1 


7(/+l) 


km 


d 


dE, 


■i n(EkmPk(wjmn  |  ek)) 


km 


J  Ekm=E^ 


therefore  (44)  can  be  rewritten  as3 

d 


dE , 


X!  ^o(w jm'n)p(I\k  I  jm'n) 


I  j,m',n 


X  ln(£te'ft(W,m'n  I  ®k))~JEkm 


=  0. 


*-=££T1> 


It  should  be  noted  that  if  GEM  is  used  to  update  Qk 
according  to  (38),  then  the  inequality  in  (47)  is  still 
true  (given  a  small  enough  stepsize)  since  (38)  implies 
gradient  ascent  of  the  function  in  curly  brackets  of 
(46).  Equation  (47)  can  be  modified  by  summing  both 
sides  over  k ,  using  (37)  to  simplify,  then  cancelling 
terms.  The  result  is 

X  Po^jmn)p(,\k  I  j^^km^Pk^jmn  I 

j,k,m,n 

>  X  Po^jmn)p(,)<<k  I  imw)ln(£'2ft(w;m„  I  e^)). 


(45) 

The  rule  governing  the  update  of  <dk  is  given  by  (36), 
which  can  be  rewritten  as 

-g|-  |  X  Po(wj,nJpU)(k  I  jm'n) 

^  I  j,m',n 


Thus,  the  sum  of  the  second  and  third  terms  on 
the  right-hand  side  of  (43)  is  >  0,  which  means 
LL(0(/+1))  —  LL(0(/))  >  0.  This  implies  convergence 
to  a  local  maximum  of  LL. 
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x  \n(Ekm,pk(Wjm,n  |  Qk))-JEkm 


=  0. 

et=o‘,+1) 


(46) 

Observe  that  (45)  and  (46)  contain  the  same 
expression  within  the  curly  brackets,  which  is  a 
function  of  Ekm  and  (~)k  for  indices  k  =  1,2 
and  m  =  1,2 Furthermore,  (45)  and  (46) 
imply  the  function  within  curly  brackets  has  an 
extremum  at  the  point  (Ekm  =  Ep^\Qk  =  0[/+1)). 

By  computing  the  second  derivatives  of  the  curly 
bracketed  function  with  respect  to  Ekm  and  Qk  (not 
shown  here),  it  is  straightforward  to  show  that  the 
point  (Ekm  =  Ek+]>,Gk  =  ©f+1))  represents  a  global 
maximum,  given  the  Gaussian  components  used  in 
our  model,  and  given  that  P1'1  has  been  fixed  using 
parameter  values  from  the  previous  iteration.  Thus, 


X  Po(^jmJPa)(k  I 


x  ln(£w1}ft(w 


jmn 


> 


X  Po(wjm'n)p<l>(k  I  jm'n) 


xln(<>,(wim,J0f))-y42 


(47) 


3 It  is  important  to  recall  that  P{,\k  \  jmn)  is  fixed  using  the 
parameter  values  and  (-)<k  >  which  were  computed  in  the 
previous  /th  iteration  (refer  to  (15)).  Therefore  P,,](k  \  jmn )  is  not 
a  function  of  the  variables  E,  and  Qk  which  we  seek  to  optimize 
in  iteration  /  +  1 . 
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